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Spectra and radial flow at RHIC with Tsallis statistics in a Blast- Wave description 
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We have implemented the Tsallis statistics in a Blast- Wave model and applied it to mid-rapidity 
transverse-momentum spectra of identified particles measured at RHIC. This new Tsallis Blast- 
Wave function fits the RHIC data very well for pr <3 GeV/c. We observed that the collective flow 
velocity starts from zero in p+p and peripheral Au+Au collisions growing to 0.470 ± 0.009(c) in 
central Au+Au collisions. The (q— 1) parameter, which characterizes the degree of non-equilibrium 
in a system, changes from 0.100 ± 0.003 in p+p to 0.015 ± 0.005 in central Au+Au collisions, 
indicating an evolution from a highly non-equilibrated system in p+p collisions toward an almost 
thermalized system in central Au+Au collisions. The temperature and collective velocity are well 
described by a quadratic dependence on (q — 1). Two sets of parameters in our Tsallis Blast- Wave 
model are required to describe the meson and baryon groups separately in p+p collisions while one 
set of parameters appears to fit all spectra in central Au+Au collisions. 



I. INTRODUCTION 

Identified particle spectra in transverse momenta pro- 
vide pillars in the discoveries in relativistic heavy ion 
collisions 

EL H !, H- In Au+Au collisions at RHIC, 
identified particle yields per unity of rapidity integrated 
over the transverse momentum range have provided in- 
formation on chemical potential and temperature at the 
chemical freeze-out in a statistical analysis p], Q. The 
transverse momentum (px ) distributions of particles 
with different masses can be described in a Boltzmann- 
Gibbs Blast- Wave (BGBW) model with a compact set 
of parameters of temperature (T), flow velocity (/3 ) and 
flow profile (p) 0, @. With the BGBW model applied 
to the data in a limited pt range, a large radial flow 
(3 ~ 0.6 is obtained in central Au+Au collisions, while 
much smaller flow is obtained in peripheral Au+Au col- 
lisioins [f| Q. The description of heavy ion collisions 
using these freeze-out conditions have been used as evi- 
dence of collective motion in relativistic heavy ion colli- 
sions 0, H, H, H, [l(j (which has further allowed extraction 
of drag and diffusion coefficient of heavy quarks in the 
expanding bulk Q3, Q3, [H 03 ) , provided necessary inter- 
play with elliptic flow effect for the observed mass order- 
ing in V2 [H , and as a necessary condition for ex plai ning 
the ridge phenomenon in many models p^ . ll6l .ri7l. ll8l . [l^ |. 

However, such BGBW descriptions have limitations. 
Modeling nuclear collisions in ideal hydrodynamics is lim- 
ited to low pt (JS 1 GeV/c) because it is generally be- 
lieved that the equilibrium description fails at high pt, 
where particle production may be dominated by non- 
equilibrium or hard processes and exhibits a character- 
istic power-law tail [2(j. The blast- wave model has a 
strong assumption on local thermal equilibrium so that a 
Boltzmann distribution can be applied [|| . This results in 



an arbitrary choice of pt range of the spectra where the 
function is able to fit the data and requires low and high 
Pt cuts [BIjQ- The BGBW model also lacks non-extensive 
quantities to describe the evolution from p+p to central 
A+A collisions. For example, one would expect that the 
energy deposited at mid-rapidity for particle production 
fluctuates significantly from event to event in p+p col- 
lisions while the current BGBW treats all p+p events 
the same with a heat bath at a fixed temperature. The 
resulting finite flow velocity (3 ~ 0.2 in p+p collisions 
from the same fitting procedure obscures the interpreta- 
tion that collective flow is large and unique in A+A colli- 
sions @. In Au+Au collisions, the fluctuations at initial 
impact due to Color-Glass Condensate (CGC) formation 
or individual nuclcon-nucleon collision may not be com- 
pletely washed out by subsequent interactions at cither 
the QGP phase or hadronic phase [2l|, [22|, [23[ . All these 
effects leave footprints in the spectra at low and interme- 
diate pt {pt ^ a few GeV/c). For example, the volution 
of <p and K* spectra from an tjit power-law (Levy) func- 
tion in p+p and peripheral Au+Au collisions to an tut 
exponential (Boltzmann) function in central Au+Au col- 
lisions can be clearly observed @, Hf| HE H3; where 
rriT = \J tuq + p T is the transverse mass of particle with 
mass nio at a given pt ■ 

With its development and success of Tsallis statistics 
in dealing with complex systems in condensed matter, 
many authors have utilized Tsallis statistics to under- 
stand the particle production in high-energy and nuclear 
physics [23, [H, Ho, El H3] ■ Although the implications 
and understanding of the consequences of such an appli- 
cation are still under investigation, the function is rela- 
tively easy to understand. The usual Boltzmann distri- 
bution in an mj- exponential form is re-written as an my 
power-law function: 
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where the left-hand side is the invariant differential par- 
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ticle yield and q is a parameter characterizing the de- 
gree of non-equilibrium. The distribution can be derived 
from the usual procedure in statistical mechanics, start- 
ing from a non-equilibrium g-entropy [H|. If particle 
production is not distributed at a fixed temperature but 
rather as a system where T varies with total energy (E) 
as T oc (q — 1)E, this will result in a negative binomial 
distribution (NBD) for particle number and temperature 
fluctuations. Their relations to q are the skewness in 
NBD k = l/(q — f) and the fluctuation of temperature 
cr 2 (l/T) = (q-l)(l/T) 2 [H. Regardless of the physical 
interpretation, Eq. [T]docs provide a necessary power-law 
behavior at high px and an exponential behavior at low 
px ■ This is exactly what has been observed in p+p 
and peripheral Au+Au collisions at RHIC for K* and cj>, 
whose measurable px range can reach low and high px 
[MES EE El- When q -> I, Eq. [I] becomes the familiar 
Boltzmann distribution again. 

In addition to the features of the spectral evolution 
necessary to describe the px spectra observed at RHIC, 
the physical interpretation of the results can provide 
quantitative insight into the non-equilibrium processes 
in relativistic heavy ion collisions and whether the sys- 
tem has thermalized in central Au+Au collisions. We 
expect that individual nuclcon-nucleon collisions inside 
a nucleus-nucleus collision are in non-equilibrium at the 
initial impact and/or have a large energy fluctuation or a 
large (q — 1) value, which produce a large power-law tail 
in the mx spectra. Alternatively, in a CGC scenario, the 
system as a whole produces a strong color field with large 
fluctuations [l?], EE El] . This can be treated as many hot 
spots in a nucleus- nucleus collision at the initial impact. 
In a viscous hydrodynamic evolution, the hot spots are 
smoothed (dissipated) into producing collective flow, cre- 



ating more particles and increasing temperature [361 ] . It 
has been argued in the Tsallis statistics that the increase 
of temperature and flow velocity during the evolution is 
connected to the decrease of (q— 1) by (shear and bulk £) 
viscosity in linear or quadratic proportion [20L E3| > such 
as T = To + a£,(q — l) 2 . This can provide quantitative 
insight into the bulk viscosity, which is predicted to peak 
at the phase transition and is much larger than the shear 
viscosity [13, El]. 

In this paper, we present the procedure of implement- 
ing Tsallis statistics in the Blast- Wave model (TBW) 
and use it to fit the identified particle spectra at mid- 
rapidity at RHIC. The physics implications are discussed 
and preparation for future work is also presented. Good 
TBW fits can also provide a practical experimental tool 
to extract particle yields (dN/dy) by extrapolating to 
unmeasured kinematic ranges since most of the experi- 
mental measurements only cover a limited px range for 
any given particle. 
II. IMPLEMENT TSALLIS STATISTICS INTO 
BLAST- WAVE MODEL 

To take into account collective flow in both longitudi- 
nal and transverse directions in relativistic heavy ion col- 
lisions, a simple Tsallis distribution needs to be embed- 
ded in the framework of hydrodynamic expansion [29| . 
We follow the recipe of the Blast- Wave model provided by 
Schnedermann et al. [||. The formula has been adopted 
by many authors to implement a fit to the data in rel- 
ativistic heavy ion collisions [S @, [IS El]- It is rela- 
tively trivial to change sources of particle emission from 
a Boltzmann distribution to a Tsallis distribution in the 
Blast- Wave model. 
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where p = tanh _1 (/3 ;! (-^)' 1 ) is the flow profile growing as 
n-th power from zero at the center of the collisions to 
(3 S at the hard-spherical edge (R) along the transverse 
radial direction (r), and (3 = (3 S /{1 + l/(n+ 1)) is the 
average flow velocity. We have used n = 1 in this study. 
However, the integrations after the replacement are hy- 
pergeometric functions and the integrals over rapidity (y) 
and azimuthal angle (0) cannot be decoupled into two 
Bessel functions as in BGBW Jf|. The computing pro- 
gram typically used to provide conventional Blast- Wave 
fits to RHIC data [11] has been modified to numerically 
calculate the integration in the above spectral function 
and fit the reslting distributions to the data. There are 
a few assumptions which do not change when a Tsallis 
distribution replaces a Boltzmann distribution: 



1. Bjorken longitudinal expansion is assumed so that 
the measured particle yield does not depend on ra- 
pidity due to the integration over rapidity of the 
source @- This is approximately true at mid- 
rapidity at RHIC or LHC [H. On the other hand, 
this assumption can be lifted in future analyses 
with a more complicated integration provided an 
emitting source along rapidity is known. 

2. Isotropic emission in azimuth is assumed for each 
local source. However, the distribution of the 
source can have an azimuthal dependence in real- 
ity [4ll ]. In the future, this can be implemented as 
has been done in the conventional BGBW model 
to fit the azimuthal dependence of HBT radius and 
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the identified particle elliptic flow 

3. The emission source everywhere has the same den- 
sity and degree of non-equilibrium (g) at the time 
of kinetic freeze-out. This may not be true since 
the high-pr particles (jet) tend to have surface 
emission [42|, |43[. This kind of corona effect has 
been implemented in many other models and can 
be adopted in this framework as well. 

4. Resonance decay contributions to the stable parti- 
cle yields have been treated as part of the source 
emission. The detailed decay kinematics and its ef- 
fect on the spectra have been studied in Ref. (til.l^. 
Incorporation of resonance effects on spectra can be 
made in future improvements. 

The goal of this paper is to provide a first implementa- 
tion of the TBW model to exercise a trial case with the 
RHIC data. Future work can change the above assump- 
tions to resemble more realistic conditions by comparing 
the assumptions to data and hydrodynamic calculations. 
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FIG. 1: (Color Online) Identified particle transverse momen- 
tum spectra in Au+Au collisions at y^s^" = 200 GeV inO- 
10% central (a) and in peripheral 60-80% collisions (b). The 
symbols represent experiment data points. The solid curves 
represent the TBW fit. 



III. RESULTS OF FIT TO RHIC DATA 

The STAR Collaboration has published a series of par- 
ticle spectra at mid-rapidity. The most complete set is 
for p+p and Au+Au collisions at ^/s^ = 200 GeV. The 
identified particle spectra include tt ± 1 Ks, K*, p, <f>, 
A, 5, p, A, and S flSS^HSHSHlii- The 
invariant differential yields were measured as 2 irpTdp T dy m 
In the new TBW model, three parameters are common 
for all particles: temperature T, non-equilibrium param- 
eter and maximum flow velocity (3 S = /3(1 + n/2) 



where n = 1 and average flow velocity j3 is bounded to 
the range [0,0.7] [3^| to aid in fit convergence and avoid 
non-physical results. An additional parameter provides 
the overall normalization of dN/dy for each species. We 
choose the Minuit in Root [5(j to perform a least-x 2 fit 
used in ref. Figure [1] shows the pt spectrum data to- 
gether with our fit results in two selected centrality bins 
(0-10% and 60-80%) in Au+Au collisions. The fit pa- 
rameters and x 2 /DoF are tabulated in Tab.|TJ As stated 
earlier in our model's third assumption, surface emis- 
sion could become important at high px; we limit our 
fits to px < 3 GeV/c to avoid this region, which still 
extends the fit range well beyond previous BGBW fits. 
The curves from our model generally describe the data 
very well, especially in central Au+Au collisions. For 
peripheral Au+Au collisions, the meson spectra are well 
described by the model while the baryons are in general 
over-predicted at higher px ■ On the other hand, the 
X 2 /DoF show good fits in all cases. The main results 
are: 

1. (q — 1), a measure of the degree of non-equilibrium, 
decreases by a factor of 5 from 0.086 to 0.018. This 
means the power in the rriT power-law increases 
from about 12 to 56, attaining an almost Boltzman 
distribution. 

2. T, the average temperature of the local source, 
shows a small increase from 114 MeV to 122 MeV. 
This trend is in contrast to the conventional BGBW 
result, where a decrease of temperature was ob- 
served [39j . 

3. /?, the average flow velocity, increases from in pe- 
ripheral to 0.47c in central Au+Au collisions. That 
the minimum x 2 /DoF is found at the lower bound 
of for peripheral collisions indicates that cither 
the model is incomplete (some approximations may 
not be sufficiently true), or that no flow has devel- 
oped in peripheral collisions within the context of 
the model. This also coincides with a large q-1, 
indicating a very non-equilibrated system if the de- 
scription as a unified system applies. 

Figure [2] shows the temperature and flow velocity ver- 
sus (q — 1) for Au+Au collisions. Each shaded region rep- 
resents a one-cr contour from the error matrix obtained 
from the TBW fit for a given centrality. The dependence 
is clearly non-linear and has a negative correlation. There 
is a jump of flow velocity from zero in p+p and 60-80% 
Au+Au centrality to 0.28 at 40-60% Au+Au centrality, 
coinciding with the transition behavior in several other 
observables [H[ . We fit the distributions with a quadrat- 
ics and obtain T = (0.123 ± 0.0014) - (1.2 ± 0.4) (g - l) 2 
and/3= (0.49 ± 0.01) - (61 ± 5)(q - l) 2 , as shown in the 
figure. 

Since the TBW model can be used to describe systems 
at non-equilibrium, it is natural to extend the fit to p+p 
collisions. However, a very poor x 2 /DoF was obtained if 
we include all of the mesons and baryons in a common 
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TABLE I: Values of parameters from TBW fit to identified 
particle transverse spectra in Au+Au collisions of different 
centralities and in p+p collisions at RHIC. Quoted errors are 
quadratical sum of statistical and uncorrelated systematic er- 
rors. The limits of /3 is set to [0,0.7]. 
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FIG. 2: (Color Online) The fit parameters T and f3 as a func- 
tion of (q — 1). Each block is one-a contour from the error 
matrix of the TBW fit for a given centrality of Au+Au colli- 
sions. 



fit. Instead, two separate groups of mesons and baryons 
show good fits. Figure [3] shows the results of the fits 
together with the data points. In both cases, the flow 
velocity was set at the lower limit of (3 = as was also 
independently verified if j3 was set to be a free parameter. 
In Fig. [3Ja, the proton and anti-proton spectra are pre- 
sented together with the predicted curves from the TBW 
fit to the meson group. The need to separate mesons and 
baryons can be seen from spectral shape of the proton 
data, which matches more closely that of the kaons than 
the K* (steeper at low px ) despite being closer in mass 
to the K* . The TBW model restricts spectral shapes to 
vary monotonically with mass for fixed (T, q, /?), so there 
is no allowance for this artifact. Weak-decay effects from 
A are investigated based on the results from the baryon- 
only fit and A— decay kinematics and do not explain this 
large steepening. On the contrary, a significantly smaller 
q flattens the spectra at low px and softens them at high 



Pt , which is necessary and sufficient to achieve a good 
fit for the baryons (p, A, S) and anti-baryons. This im- 
plies that not only mass plays an important factor in the 
particle yield at a given px , the particle species also 
significantly affects the outcome in p+p collisions. This 
characteristic baryon versus meson grouping in p+p data 
has been seen previously in mx scaling analyses of the 
same data 44|, and our results confirm this observation 
that baryon number plays an important role in hadron 
production in p+p collisions. 



IV. DISCUSSIONS AND OUTLOOKS 

Modifying the Blast-Wave model to utilize Tsallis 
statistics instead of the conventional Boltzmann-Gibbs 
statistics has allowed high quality fits (see x 2 /DoF in Ta- 
ble|T| over a broader transverse momentum range and has 
altered the conclusions which can be drawn from the fits. 
The extended px range is enabled by the Tsallis statistics' 
capacity to evolve from an mx exponential source into a 
power law through increasing q values, though the phys- 
ical interpretation of this statistical model in the context 
of high energy nuclear collisions remains to be fully un- 
derstood. 

For central Au+Au collisions at ^/s NN = 200 GeV, the 
TBW fits find a value of q approaching unity, which im- 
plies results very similar to those from BGBW with a 
large radial flow velocity Q. The progression of the 
TBW fit parameters with centrality is generally smooth 
but is strikingly different from BGBW fits, culminating in 
a notably non-equilibrium description of p+p and periph- 
eral Au+Au collisions where a preference is found for zero 
collective flow. Qualitatively, increasing centrality pro- 
duces a strong increase in flow velocity, a mild increase 
in temperature, and a dramatic decrease in (q— 1). This is 
consistent with a picture of increased thermalization with 
centrality but disfavors complete thermal equilibrium in 
all systems, a requirement for applicability of the BGBW 
model. The BGBW model appears to translate the non- 
equilibrium features, indicated by our model, into non- 
zero collective flow velocity and higher temperatures in 
p+p and peripheral Au+Au collisions 0, @] 

Hydrodynamics with space-time evolution from an ini- 
tial condition [36| is so far the most realistic simulation 
of what happens in relativistic heavy ion collisions. How- 
ever, even with implementation of initial conditions and 
an interface to hadronic cascade models at late stage of 
the evolution, it is hard to obtain an intuitive picture, in 
contrast to that offered by an analytical parametrization 
in the Blast- Wave model. It seems unlikely that hydro- 
dynamics is applicable at all for p+p and very peripheral 
A+A collisions at RHIC. Being able to provide a system- 
atic comparison between p+p and central A+A collisions 
in one model framework is still valuable and in some cases 
may be necessary. 

We have additionally observed that p+p spectra con- 
tinue to be well described by our TBW curves from 3 to 
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10 GeV/c, well beyond the range used for the fits. This 
implies that the origin of the power law behavior may be 
the same from low pt to high px, whether their under- 
lying connection is fragmentation, parton evolution, or 
parton interaction cross-section. However, the high pt 
behavior in central Au+Au collisions is notably different: 
the experimental spectra have a reduced power-law tail 
relative to the binary scaled p+p spectra, but still signifi- 
cantly above the extension of the fit from px < 3 GeV/c. 
In the near future, we plan to implement a corona-like 
radius-dependent (q— 1) which could accommodate larger 
power-law tails in central Au + Au collisions and allow the 
fit range to be extended to higher px ■ This could provide 
additional information on jet quenching as characterized 
by surface emission. 
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FIG. 3: (Color Online) Transverse momentum spectra of iden- 
tified mesons (a) and baryons (b) produced in p+p collisions 
at ^fs = 200 GeV. The symbols represent the experiment 
data. The solid curves represent the TBW fit. The proton 
and anti-proton transverse momentum spectra and predicted 
curves from meson fit parameters are also plotted on panel 
(a) for comparison. 

It has been argued [2(| Hi| that the dependence of 
T and (3 on (q — 1) is related to bulk viscosity. If this 
viscosity is very large at phase transition, a systematic 
study of TBW model fits from AGS to SPS and RHIC 
energies may help locate the critical point. This, how- 
ever, requires measurements of many identified particles 
to at least intermediate pr (3 GeV/c) at all energies in 
a systematic fashion. These data can be provided by 
a RHIC beam energy scan. At LHC, one expects even 
larger (q — 1) value in p+p collisions than that at RHIC 
due to increased relative contributions of hard and semi- 



hard processes; the conventional BGBW is not expected 
to be very meaningful for p+p and peripheral Pb+Pb 
collisions. If (g — 1) is not larger for LHC p+p collisions 
and a non-zero flow velocities is observed, the question 
of thermalization or even QGP in such collisions might 
be raised. The results of TBW fits at LHC will certainly 
be informative. 

V. CONCLUSIONS 

In summary, we have implemented the Tsallis statistics 
in a Blast- Wave model and applied it to sets of identi- 
fied particle spectra versus transverse momenta at mid- 
rapidity at RHIC. This new TBW function fits the RHIC 
data quite well for px <3 GeV/c. We observe that the 
collective flow velocity starts from zero in p+p and pe- 
ripheral Au+Au collisions and rises to 0.470±0.009 c in 
central Au+Au collisions. The parameter [q — 1), which 
characterizes the degree of non-equilibrium in a system, 
changes from 0.100 ± 0.003 to 0.015 ± 0.005 systemati- 
cally from p+p to central Au+Au collisions, indicating 
an evolution from a highly non-equilibrated system in 
p+p collisions toward an almost thermalized system in 
central Au+Au. TBW fits using all species from mul- 
tiple ccntralitics demonstrate a quadratic dependence of 
the temperature and collective flow velocity on (q — 1), 
diverging into two different fits for mesons and baryons 
necessary to best describe p+p spectra. 
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